Modelling the effect of gap junctions on tissue-level cardiac 

electrophysiology 



Doug Bruce 



Pras Pathmanathan 



Jonathan P. Whiteley 



douglas .bruceOoriel . ox . ac .uk pras . pathmanathanScs . ox. ac .uk Jonathan. whiteleyOcs . ox . ac .uk 
Computational Biology Group, Department of Computer Science, University of Oxford 

When modelling tissue-level cardiac electrophysiology, continuum approximations to the discrete 
cell-level equations are used to maintain computational tractability. One of the most commonly used 
models is represented by the bidomain equations, the derivation of which relies on a homogenisa- 
tion technique to construct a suitable approximation to the discrete model. This derivation does not 
explicitly account for the presence of gap junctions connecting one cell to another. It has been seen 
experimentally [Rohr, Cardiovasc. Res. 2004] that these gap junctions have a marked effect on the 
propagation of the action potential, specifically as the upstroke of the wave passes through the gap 
junction. 

In this paper we explicitly include gap junctions in a both a 2D discrete model of cardiac elec- 
trophysiology, and the corresponding continuum model, on a simplified cell geometry. Using these 
models we compare the results of simulations using both continuum and discrete systems. We see 
that the form of the action potential as it passes through gap junctions cannot be replicated using a 
continuum model, and that the underlying propagation speed of the action potential ceases to match 
up between models when gap junctions are introduced. In addition, the results of the discrete simu- 
lations match the characteristics of those shown in Rohr 2004. From this, we suggest that a hybrid 
model — a discrete system following the upstroke of the action potential, and a continuum system 
elsewhere — may give a more accurate description of cardiac electrophysiology. 

1 Introduction 

Many phenomena in biology are discrete; for example, biological tissue consists of discrete cells within 
extracellular material. When modelling a particular phenomenon or feature, different sets of equations 
can apply in the intra- and extracellular regions. In the case of modelling at the tissue or organ level it is 
impractical to model each individual cell. As a result, multi-scale techniques which consider the average 
behaviour of the problem — for example homogenisation — are used so that we may include cell-level 
phenomena into a tissue-level model whilst retaining computational tractability. 

In the case of cardiac electrophysiology, whilst computing a numerical approximation to the solution 
of the governing equations at the level of individual cells should give accurate and realistic results, it 
is computationally unfeasible even for very small regions of tissue. It is much more efficient to solve 
these models using a continuum approximation to the equations, and indeed over the last 20 years many 
have done this [8 ]. This continuum model relies on a homogenisation technique to construct a suitable 
approximation to the discrete model. The homogenisation process uses a multiple-scales method which 
uses the fact that the problem is naturally defined on two different scales. 

We now discuss relevant issues from cardiac electrophysiology and homogenisation. 
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1.1 Cardiac tissue modelling 




Figure 1: Cardiac cell structure: Guy ton and Hall, 1996, Fig. 9-2, p. 108. 

An example of the histology of cardiac cells is shown in Figure [T] These myocardial cells are roughly 
cylindrical and are packed together in an irregular three-dimensional pattern. Each cell is electrically 
connected to its neighbours via gap junctions, which are small channels through which ions may flow. 
The cell structure is surrounded by an extracellular matrix (ECM) and the two regions are separated 
by the cell membrane. In the discrete formulation of the governing equations, as will be described in 
more detail shortly, Laplace's equation holds for the potentials in both the intracellular and extracellular 
spaces, with the electrical properties of the membrane taken into account in the boundary conditions. 
The solution of these equations is a substantial computational task even for a small region of tissue due 
to the fine-scale structural detail of cells. As a result, it is computationally desirable to homogenise the 
microstructure to allow us to pose the problem as a continuum rather than a discrete problem. This led 
to the proposal of the bidomain model by Tung in 1978 [12j in which the homogenised potentials are 
solved for via two reaction-diffusion equations using averaged conductivities, where the precise form 
of averaging arises from the homogenisation technique. A sketch of the derivation of the bidomain 



equations is presented later in this document (see Section 2.2 ), and as discussed earlier it is the generally 



accepted model of cardiac tissue behaviour. A more formal derivation of the model can be found in (5J. 

Whilst realistic simulations using the bidomain model have been demonstrated, these studies tend to take 
the bidomain equations as being physiologically correct in all circumstances, and empirical parameters 
are measured by matching the results of the continuum simulations with experimental data. However, 
there has been little rigorous testing of the validity of the derivation of the bidomain equations; in partic- 
ular, concerning the homogenisation technique used to average microstructural quantities. 

In addition, the effect of gap junctions on propagation is often ignored during the homogenisation pro- 
cess. As previously mentioned, gap junctions are the means by which cells are electrically connected to 
one another. Whilst they allow the signal to be conducted, they do so with more resistance than is given 
by the interior of the cells. When considering action potential propagation at cell-level, we therefore 
expect to see the conduction velocity reduce as the wave passes through the gap junction. Indeed, this is 



supported by experimental data as shown in [ 1 1 ] — the spatial form of the action potential is 'stepped', 
with the reduced conductivity in the gap junction causing a steep jump in membrane potential between 
one side of the gap junction and the other. The repercussions of this for the continuum model, and how 
it affects the homogenisation process, will be discussed later. 
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1.1.1 Hybrid models of cardiac electrophysiology 

In recent years, models based on hybrid automata have been developed to model networks of excitable 



cells, with an initial focus on the temporal morphology of the action potential |13|. The model was 
further refined to study more specific conditions associated with action potential morphology such as 
early depolarisation p4| , spiral waves (3j, and tachycardia [4]. Whilst these models allow efficient and 
precise analysis of the conditions required for such phenomena, they are not focused on the incorporation 
of the cell microstructure and the corresponding effect on action potential propagation. 

1.2 Homogenisation 

Whilst the initial formulation of the bidomain model stated the equations in the same form as they 
remain today (as laid out in [ 12]), no rigorous derivation was given at that time that followed from basic 
physical principles and the cell-level properties of the tissue — the macroscopic bidomain equations were 
not directly connected to the microstructure of the tissue. The advent of more formal homogenisation 
methods (2j paved the way for a first attempt at a rigorous derivation of the bidomain equations by 
Neu & Krassowska in 1993 |7]. They used a multiple-scales asymptotic expansion technique to convert 
the microscopic problem, formulated in terms of the pointwise potentials, to a macroscopic problem 
formulated in terms of the leading order potential averages. 

This technique was then modified by Keener & Panfiiov in 1996 (5J, who corrected a couple of errors 
in the homogenisation technique used. This boiled down to generalising the work of Neu & Krassowska 
using more realistic tissue geometries. The error corrections played no part during normal action po- 
tential situations, but had a significant effects when examining the result of the application of a large 
current stimulus such as is seen during defibrillatory shocks, with the result being that the mechanism for 
defibrillation was clarified. It is a version of this derivation that will be presented later in this document. 



Following this, a recent paper by Richardson & Chapman [10] augmented the above derivation by in- 
troducing a co-ordinate transformation that allows for variable tissue structure and ultimately a set of 
bidomain equations in which the conductivity tensors systematically account for deformation of the tis- 
sue and the orientation of the cells. 

1.3 Analysis of the homogenisation process 

Despite all this work having been done on the derivation of the bidomain equations, looking at multiple- 
scales methods and complex geometry, one unproven assumption remains. In the derivations, the macro- 
and microscales of the problem are related by a dimensionless parameter, usually denoted e, that is equal 
to the ratio of the typical lengthscale of a single cell to the lengthscale over which the solution varies. 
In fact, the denominator of e is stated as the 'natural' lengthscale of the fibres in pj, which has been 
translated in pOj as the 'typical lengthscale of the cardiac tissue', but as we have stated above it is more 
correct to translate it as the typical solution lengthscale, so that e is equal to the ratio between microscale 
and macroscale coordinates. This parameter is assumed to be small, and much of the homogenisation 
process involves taking the limit e — > 0. 

However, it is clear that in the case of a steep propagating wavefront, this will not be the case. As pre- 
viously discussed, when the action potential passes through a gap junction the steepness of the upstroke 
portion of the wave will be significantly increased. We therefore question whether, in such circumstances, 
the homogenisation process used in the derivation of the bidomain equations remains valid. 

Furthermore, it will become apparent when presented with an explicit version of the homogenisation 
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process in Section 2.2 that the 'stepped' action potential seen in the presence of gap junctions cannot 
be captured using such a continuum model — not only does the homogenisation process enforce that 
the macroscale conductivity tensors are independent of space at a cellular level, but the fundamental 
principle of the continuum model is that the fine-scale structure of cells, including gap junctions, are 
only accounted for by their average effect on the system. 

1.4 Aims and Outline 

We have described how the presence of gap junctions in cardiac tissue has a major effect on the propa- 
gation of the action potential and the associated conduction velocity. We therefore wish to consider the 
most accurate method to incorporate such structures into a simplified discrete model of cardiac electro- 
physiology. 

We begin in Section [2] by describing the conventional discrete and continuum models of cardiac elec- 
trophysiology. Then, in Section [3} we consider the physiological properties of gap junctions and relate 
these to the models, setting up a more detailed discrete model that accounts for gap junctions, as well 
as deriving the corresponding continuum equations, which we will refer to as the modified bidomain 
equations. We then compare the results of these continuum and discrete models in Section[4j asking how 
well the solutions compare both to each other and to those seen experimentally. We will observe that 
whilst the bidomain equations adequately represent the behaviour of a discrete model that neglects gap 
junctions, when gap junctions are included the continuum model wavespeed begins to fail to match the 
discrete model wavespeed, and that the discrete model better captures action potential profiles observed 
experimentally. 

2 Modelling cardiac electrophysiology 

The derivation of the models below applies to general periodic structures. However, in this paper we are 
restricting our considerations to 2D models, and so we will now present a simplified representation of 
the cell structure in two dimensions. This will be the geometry upon which our analysis and simulations 
will take place. 

To arrive at our simplified representation of the cells, we assume that the intracellular space in each peri- 
odic subunit is rectangular in shape, with no coupling of the cells in the off-fibre direction, thus restricting 
our attention to propagation in one spatial dimension. This leads us to the schematic representation of 



Figure 2a 



This representation consists of a sheet of cells that are connected into fibres along the x-direction, with 
each fibre separated in the y-direction by extracellular matrix. We create a periodic subunit of the domain, 
labelled H, which contains both intracellular and extracellular portions, and this can be seen in more 



detail in Figure 2b Here, we are assuming that all cells have the same dimensions. 



2.1 The discrete model 



In this section, we are referring to the global problem as given in Figure 2a 
2.1.1 Intracellular space 

In the intracellular space £2,-, the current is given by i c = — a, V(j>i, where 0,- is the intracellular potential and 
<7i is the (scalar) conductivity of the intracellular space. The form of the current comes from Ohm's law 
which states that current is the product of conductivity and the gradient of the potential [6]. Conservation 
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(a) A representation of cardiac cells in 2D. (7,- and 
O e are, respectively, the intracellular and extracellular 
conductivities, with 0,- and <j> e the potentials. The re- 
gion fl represents a periodic subunit containing both 
intracellular and extracellular space. 




L 

(b) A single periodic subunit of cell and ECM, as repre- 
sented by fl in Figure [2a] L represents the length of our 
subunit, with h its height and hi the height of the intra- 
cellular portion. The intracellular domain is labelled £2,-, 
and the extracellular domain Q. e . The membrane between 
the two is labelled d£l m . 



Figure 2: The 2D representation of cardiac cells 



of current in the intracellular space therefore gives us 

v-(ojVfc) = o, xea,. (i) 

The boundary condition representing flux of current is then 

-OjVft • n = I m {x), x e da m , (2) 

where n is the outward pointing normal, that is, the normal pointing from the intracellular space into the 
extracellular space, I m is the transmembrane current, i.e. the current flowing from the intracellular space 
into the extracellular space, and dQ. m is the boundary between the intracellular and extracellular spaces. 

2.1.2 Extracellular space 

Similarly, in the extracellular space Q. e we have 

V-{a e V<j) e ) = 0, xeO, (3) 

where <p e and o e are the extracellular potential and conductivity respectively. The transmembrane current 
I m will now flow into the extracellular space, and so the boundary condition here becomes 

a e V(j) e ■ n = /,„ (x) , x e da m , (4) 

where dQ. m is the same boundary and n the same normal as in ([2]), so that the normal still points from 
the intracellular to the extracellular space. 

2.1.3 The transmembrane current 

For the transmembrane potential, defined on dQ. m by v m = (pi — ty e , the transmembrane current is given 
by 

Im — C m — h /('on(v m , u) , (5) 
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where C m denotes the membrane capacitance and u(x,f) consists of various ionic concentrations and 
gating variables, which are determined from a cell model usually represented by a system of ordinary 
differential equations (ODEs). This comes from modelling the cell membrane as a capacitor, and as the 
capacitance of an insulator is defined to be the ratio between charge and potential we have C m = Q/v m . 

Since the current is given by the rate of change of charge, i.e. dQ/dt, it follows that the capacitive current 
is equal to C,„^f-, assuming that C m is a constant property of the material. In addition to this capacitive 
current there will be an ionic current created by the flow of ions through the membrane, and this is 
denoted The total transmembrane current will be the sum of the capacitive and ionic currents, and 
is thus given by ([5]). 

As stated previously, the quantity Z ;on represents the sum of the currents formed by the flow of charged 
ions across the cell membrane. These currents are due to the membrane possessing pore-forming pro- 
teins, known as channels, which allow the passage of specific ions, for example sodium (Na + ) and 
potassium (K + ), down their electrochemical gradient. The specifics of the form of the current caused by 
these channels can be found in [6j. It is assumed that / i0 „ is defined in the same fashion at all points on 
the cell membrane, and that there are a sufficiently large number of ion channels along the cell membrane 
for us to model l\ on as a continuous function. 

To summarise, in a discrete framework we will be solving the equations [([T]>,([3]),(|5])] for the unknowns 0,- 
and (p e , subject to the boundary conditions [pKQ], as well as specifying that solutions and current are 
continuous between one cell and the next. 



2.2 The continuum approximation 

To model the cardiac tissue as a continuum, we use a homogenisation technique based on the assumption 
that the lengthscale of the solution to the governing equations is much larger than the length of an 
individual cell. We therefore define a "fast" variable 

1 

z = -x, 
£ 

where 

length of a single cell 



lengthscale of the solution 
and it is thus assumed that £<1. Considering the intracellular space to begin with, we can now write 
as a function of both x and z and seek a solution by expanding in powers of £, so that 

fax,z) =< £,'(x) + £fa (x,z)+£ 2 /2 (x, Z) + ..., (6) 

where fa , fa, . . . , are periodic in z with zero mean. A full version of the derivation can be found in (71 or 
(5j, but to summarise, by substituting (|6]) into the discrete governing equations and boundary conditions, 
equating powers of £ and using the periodicity of fa and fa, the equation for the intracellular potential 
can be written 

V x - (EfV,*,-) = -!- / I m dS z , (7) 

Vcell J dQ.„, 

where V C eii is the volume of our periodic subunit and is the macroscale intracellular conductivity tensor 
of the problem, given by 

1 Ufz+^W (8) 



Vcell JCli V dz 



D.A. Bruce, P. Pathmanathan & J.P Whiteley 



7 



The functions Wj are periodic in z with zero mean, and satisfy 

V z -(a,-V z Wj) = -^, xefli, V x Wj-n=-nj, x€da m , 7 = 1,2 

The same method, applied this time to the extracellular space, will give us 



where analogously 



V x • (I e V x <D e ) = - J- / l m dS z , (9) 
Vcell JdSlm 



1 /■ / dW e 
E * = i7- / M'+^T" ) d ^' ( 10 ) 



The functions now satisfy 

V z -(a e V z W/) = -|^, xeil e , y z wj-n = nj , xedn m , j = 1,2 

and are periodic in z with zero mean. Thus, the tissue-level conductivity tensors and £ e will depend 
on the domain shapes £2, and Q e , the volumne of our periodic subunit V ce u, the micro-level conductivity 
scalars a, and o e , along with any quantity that will change the functions W^jv 

2.2.1 The bidomain equations 

If we now write the transmembrane potential v m as a power expansion in e, so that 

V,„ = V m (x) + £V,„i (X, Z) + £ 2 V m2 (X, Z) + . . . , 

we may express ([7]) as 

V x • (E,V,4»i) = ^- / C m %+ 7i OT (v m ,0 dS z , 

= 2 ( c m^r + 7T- / /,- on (V m (x) + £v m i(x,z)+e 2 v m2 (x,z) + ...,0 d5 z ) , 

where 5 m is the membrane surface area and % = Sm/Vcell- If we make the assumption that we may ignore 
the contribution from e 1 and higher order terms, i.e. that 

y j hon {V m (x) + ev m i (x, z) + e 2 v m2 (x,z) + ...,t) dS z = I ion (V m (x), t) 

"m J dQ. m 

which also physiologically implies that we are taking there to be sufficiently many ion channels in each 
cell to model Ij on and a continuous function as mentioned previously, we can write 

V x • (EfV,*,) = X ( C '"^ +hon{V m ,t)\ , (11) 



and similarly (|9]> becomes 



V x .(V.V x $ f) - _ z ( Cm ^+I ion (V m ,t) ) . (12) 



We now use the fact that 0, = V m + e to eliminate </>, from ([TT]i, and denote the transmembrane potential 
V m simply as V to avoid any future subscript confusion, to write ( [TTj ) and ( [12] ) in the more familiar form 
of 

dV 

XQ n -^ = V x • (E,-V«(V + fc)) - ^ on , (13a) 
V x .((Z, + E B )V x ^ + EfV,V)=0. (13b) 
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Appropriate boundary conditions, imposing zero flux on the boundary of the entire domain, are 

-E,-V(V + ^)-n = 0, L e V<p e -n = 0. 

These are what we call the bidomain equations, in the absence of any external stimuli. 

3 Explicitly incorporating gap junctions into the models 

3.1 Physiology of gap junctions 



Closed Open 




Figure 3: A representation of gap junctions, Mariana Ruiz, 2006. 

As seen in Figure [3} the gap junctions form a channel that directly connect two adjacent cells, allowing 
molecules and ions to pass through it. These junctions are abundant in cardiac muscle, and allow direct 
electrical signalling from one cell to the next. When looking at an entire cell, we can approximate 
the collection of individual gap junctions by one continuous domain whose conductivity o g takes into 
account the density of gap junctions present, the average fraction that are open or closed during signal 
propagation and the underlying conductivity of the material. 

We then notice that, on the boundary between such a gap junction domain and the intracellular spaces 
that it connects, both ions and electrical signals are free to flow, and so we will have continuity of flux 
across the boundary when considering electrical potential. This allows us to treat the gap junction as 
part of the intracellular space with a different conductivity, so that the intracellular conductivity becomes 
a function of space, which naturally corresponds to imposing continuity of potential and flux on the 
interface if the gap junction had been treated as a separate compartment. 

On the interface between gap junction and extracellular space, we note that there are likely to be dif- 
ferences between the properties of this interface and that of the cell membrane. Unlike on the cell 
membrane, a gap junction will not allow ion transport between itself and the extracellular space. In ad- 
dition, the capacitive properties of the material will be different from those of the cell membrane. These 
differences will have an impact on the form of the transmembrane current I m that will be defined on the 
interface between gap junction and extracellular space. 
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Figure 4: Our periodic subunit (see Figure 2b) modified to include a gap junction. It is modelled as a 
region of width 8 at one end of the cell with different conductivity ( G g ) to that of the cytoplasm ( d). 



From our discussion above, we adapt our cell geometry to include gap junctions as shown in Figure |4] 
We include a thin region of width 8 at one end of the intracellular space in which the conductivity is 
given by G g . The overall intracellular conductivity is now denoted by a = cj(x), where 

I Gi x£ cell, 

o(x) = { '. . (14) 

I G g x E gap junction. 

In addition, we adapt the form of the transmembrane current on the boundary as follows: 

Im = S c "f +Ii0 " xGce11 ' (15) 
{ c s Tt + h^on x E gap junction. 

where I g is a boolean switch that turns ionic flow across the membrane on or off, and c g is the capacitance 
of the gap junction membrane. We may then alter the properties of the gap junction membrane by 
specifying c g according to one of three cases: 

• c m , treating it as if it were the same material as the cell membrane 

• c g / c m , treating it as a capacitive material with its own properties 

• 0, treating it as a fully insulative material. 

Whilst we expect that it is correct to take I g = and c g ^ c m in the above system, we will leave both 
parameters undetermined in the forms stated so that we can explore the effect that they have on results 
of simulations. 

To summarise, the discrete equations will be modified as such: the intracellular conductivity will now 
take different values in the cell and the gap junction, and the transmembrane current will have a different 
formulation on the cell membrane and the gap junction membrane. 

3.2.1 A modified formulation of the bidomain equations 

Returning to ((7]), the right-hand term of the equation given by ^— J dilm I m dS will be split into integrals 
over the cell membrane 512, and gap junction membrane d£l g separately before being evaluated, so 
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1 



V ce ii JdQ., 



I m dS 



1 



Vcell 
1 

Vcell 



I m dS+ / I m dS 
dV dV 

CjSj— \-CvSc,^: h / IiondS-\-Io i L, on 

at Jdcii Jdci, 



dt 



f IiondS + Ig f 



liondS 



dV 



{ciXi + CgXg) + iXi + hXg)h 



(16) 
(17) 
(18) 



where c g and c, are the capacitances of the gap junction and cell membranes, S g and 5, are the surface 

s 

areas of dQ.f, and d£lj respectively, with Xg = and %i = y • O n our simplified geometry as given in 

^and^f. 



Figure [4] we will have that %i 



The final form of the modified bidomain equations will therefore be 

dV 

(ciXi + CgXg) = v x • (r ; -v x (y +&))- + hx g )h 

V x -((E l - + E e )V x ^+EiV x y)=0. 



(19a) 
(19b) 



3.2.2 Effect on the intracellular conductivity tensor 

On a generalised geometry, we see that the tensor Z, will be changed via the solutions to Q, as the 
conductivity scalar there denoted a, is modified to be a function of space. 

On our simplified geometry, we are able to write down analytic solutions for Z, in both cases. In the 
absence of gap junctions, the functions W[ and W{ given in ((9]) satisfy 

V 2 W!=0, 7 = 1,2. 
The normal is given by (0, 1) and so the boundary conditions become 



dW{ 



0, and 



~dy~ 



1, 



the solutions to which are 



W[=A h wi = -y+A 2 , 

where A \ and A2 and constants (note that they are not functions of x as W[ and W\ are periodic in x) 
Substituting this into ([8]) gives 

Vintra ( , 



Vcell V J ' 

where Vi ntra is the volume of the intracellular space. 

In the presence of gap junctions, the equations for remain unchanged, and now W[ will satisfy 
Laplace's equation in both the intracellular space and the gap junction separately with the same boundary 
condition as before. Therefore we have 



dW[ jfii 0<x<<5, 
dx |fi 2 8<x<L. 

To find B\ and B2 we integrate the governing equation for W{, given in ([9]), across x = 8 to give 

dwr° 



(21) 
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and therefore that 

a, #2 — GgBi = Gg ~ °i- 
We then use the fact that W[ has zero mean in x to give us 

dWi 



f dW 

JSli ox 



Solving for B\ and B2 and substituting into ([8]) gives us that 

v (l,l) _ Vjntra OjGgL 
2-11 — X 



Vceii 5a,- + (L-5)CTg" 

Note that, when o s = a,, this reduces to the corresponding entry in (20 1 as anticipated. It is also worth 
pointing out that, both with and without gap junctions included, we will have 

'extra I ®e 



where V ex tra is the volume of the extracellular portion of our periodic subunit. 



4 Results of simulations 



In order to investigate the effect of gap junctions on results of simulations, we solve both the continuum 
system, given by the modified bidomain equations ( fT9] ), and the discrete system, given in [QjQj©;©] 



with the transmembrane current given by ( |T5j ), using a Beeler-Reuter model for the ionic current [ 1 ] in 
both cases. We use a finite element method |9j with 320 nodes per cell in the discrete case and 80 nodes 
covering the corresponding area in the continuum case, and use the PETSc library 
(http : / /www . mcs . anl . gov/petsc) to solve the resulting linear systems. 

We take an individual subunit to be of size 0. 1 mm by 0.02 mm, with the intracellular portion 0. 1 mm by 



0.01 mm (so that, in the notation of Figure 2b we have L = 0.1 mm, h = 0.02 mm and h\ = 0.01 mm), 
using a 100 cell by 2 cell region. We simulate 50 ms of electrical activity, beginning the simulations 
in equilibrium so that 0,- = V eq and e = everywhere. Conduction coefficients are set at a,- = 0.175 
/iS/mm G e = 0.7 /iS/mm, and the membrane capacitance c m = 0.01 /iF/mm 2 , with these parameters 
taken from [8]. We then apply an appropriate current stimulus, dependent on our solution parameters o g , 
c g and I g , between 5ms and 10 ms to both cells on the v-axis. 

The table below summarises the different parameter sets used in each of our simulations, along with a 
verbal characterisation of what we are simulating. 



Parameter Values used in simulations 



Model 




c g 


h 


Characterisation 


Base 


0.175 


0.01 


1 


No gap junctions, models reduce to original forms 


1 


0.00175 


0.01 


1 


Gap junctions, simply reducing conductivity 


2 


0.00175 


0.01 





Gap junctions, not allowing ion transport across membrane 


3 


0.00175 


0.001 


1 


Gap junctions, reduced capacitance 


4 


0.00175 


0.001 





Gap junctions, reduced capacitance, no ion transport 


5 


0.00175 





1 


Gap junctions, fully insulating 


6 


0.00175 








Gap junctions, fully insulting, no ion transport 
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Base model: no gap junctions 



H -60 




Distance along fibre (mm) 



Distance along fibre (mm) 



(a) 



(b) 



Figure 5: Results of simulations at 15 ms and 30 msfor our base case (left), in which gap junctions were 
not modelled, and our most simple implementation of gap junctions (right) in which they were treated as 
a region of reduced conductivity with identical properties to that of the cell. 



Continuum Models (all) 



Discrete Models (all) 




Distance along fibre (mm) 



Distance along fibre (mm) 



(a) 



(b) 



Figure 6: Results of simulations at 30 ms for all models, comparing the continuum solutions (left) and 
discrete solutions (right) against one another, to see how our choice of model affects propagation. 




Figure 7: A magnification of the plots shown in Figure\6\ 
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4.1 Effect of introducing gap junctions 

In Figure [5] we take spatial snapshots of the results of our simulations at two separate time points — 15 
ms, seen towards the left of each subplot, and 30 ms, seen at the right of each subplot — for two of the 
models mentioned above, in both the discrete and continuum formulations of the problem. We compare 
our Base model in which gap junctions are not modelled, to Model 1 in which the gap junctions are 
modelled as a region of reduced conductivity whose membrane properties are identical to those of the 
remainder of the cell. 

As expected, in the absence of gap junctions the discrete and continuum models give near-identical so- 



lutions, as seen in Figure 5a In such a situation the problems highlighted in the introduction concerning 
the derivation of the continuum model are not applicable, and thus the continuum system provides an 
accurate representation of the discrete problem. However, when gap junctions are introduced it is seen 



in Figure 5b that the propagation speed in the continuum model does not match that of the discrete 
model. The conduction velocity of the wave for the discrete problem is noticeably smaller than that for 
the continuum problem, as observed by the small discrepancy in the solutions after 15ms and the larger 
discrepancy after 30ms. 

This discrepancy occurs because of the rapid spatial variation in the membrane potential in the discrete 
case as the action potential propagates through a gap junction. Here, the solution lengthscale is of the 
same order as an individual cell, and thus our key assumption when deriving the continuum model, 
that we may ignore effects at cell-level and below, is no longer true. It is therefore the case that the 
bidomain equations, when derived using the inherent cell-level parameters of our system, cannot be used 
as an accurate representation of the propagation of the action potential if we include the effects of gap 
junctions in our discrete model. 

In addition, we also notice that the discrete simulations in the presence of gap junctions display the form 



of 'stepped' action potential that is seen experimentally in [ 1 1 ]. As this is not seen in the absence of gap 
junctions, we conclude that there should be some representation of a gap junction structure in a model of 
cardiac electrophysiology in order to capture this more detailed form of the propagated action potential. 
It is worth noting that the continuum model, in fact any continuum model, is unable to replicate such 
behaviour — by its nature it cannot have quantities, in this case the intracellular conductivity, that vary 
on the level of single cells. 

4.2 Comparing implementations of gap junctions 

Having seen that gap junctions change the results of both discrete and continuum simulations of cardiac 
electrophysiology, ultimately causing the solutions of the two types of model to diverge, we wish to 
see if the precise nature of the implementation of gap junctions further affects the characteristics of 
solutions. To that end, in Figure|6]we plot the results of simulations of both continuum (left-hand figure) 
and discrete (right-hand figure) versions of each model specified in Table [4] taken at a time of 30 ms. A 
magnification of Figure [6]is given in Figure [7] 

Whilst we see a difference in the position of the propagating wave — and thus the underlying wavespeed 
— between each of the implementations of gap junctions, this change is much smaller than the initial 
change brought by the introduction of gap junctions over our Base model. This suggests that the major 
cause of the discrepancy between continuum and discrete solutions is the sharp change in conductivity 
that we have between the cell and the gap junction. 

Considering the plots in more detail, we see in Figure [JJ that switching the ionic current off on the gap 
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junction membrane — going from Model 1 to 2, Model 3 to 4 or Model 5 to 6 — slows down the 
propagated wave as expected, though by an equal amount in the continuum and discrete cases. Reducing 
the capacitance of the gap junction membrane — Model 1 to 3 and Model 2 to 4 — again slows down 
propagation, this time by a larger amount. However, further reducing the capacitance to zero — Model 
3 to 5 and Model 4 to 6 — has a negligible effect on solutions. More importantly, we can see that such 
changes in the results of the discrete model are mirrored in the continuum formulation of the problem, 
specifically the associated change in propagation speed. 

It is also clear from Figure[7]how the steepness of the propagating wave varies from continuum to discrete 
models — in the continuum case the wave is moderately steep for the entirety of the upstroke, whereas 
in the discrete case the wave is fairly shallow as it passes through each cell, and extremely steep inside 
the gap junction. This reinforces our statement that gap junctions cause rapid spatial variation in the 
potential. 

5 Conclusion 

The implementation of gap junctions into a standard model of cardiac electrophysiology causes a dis- 
crepancy to occur between results of simulations of discrete and continuum versions of the system. This 
is due to the rapid spatial variation in the membrane potential inside a gap junction that is caused by the 
concomitant large decrease in conductivity in such a region. Given this, it is not possible to model the 
contribution and effect of gap junctions on cardiac electrophysiology using a continuum system, and we 
suggest that a hybrid method — using the discrete model around the upstroke of the propagated wave, 
and a continuum model elsewhere — may enable us to retain accuracy and characterisation of solutions 
whilst increasing computational tractability. 

With regard to the precise implementation of gap junctions, we have seen that solutions will depend on 
the value of the capacitance of the gap junction membrane, with an adapted version of the continuum 
system matching the changes predicted by the discrete system. These results suggest it is important to 
have an accurate value for the capacitance of the gap junction membrane when conducting simulations 
of cardiac electrophysiology. 
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